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We formulate a conserving gapless mean-field theory for Bose-Einstein condensates on the 
basis of a Luttinger-Ward thermodynamic functional. It is applied to a weakly interacting 
uniform gas with density n and s-wave scattering length a to clarify its fundamental thermo- 
dynamic properties. It is found that the condensation here occurs as a first-order transition. 
The shift of the transition temperature ATc from the ideal-gas result To is positive and given 
to the leading order by ATc = 2.'i'ian^^^To, in agreement with a couple of previous estimates. 
The theory is expected to form a new theoretical basis for trapped Bose-Einstein condensates 
at finite temperatures. 
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Much effort has been devoted recently to calculating 
the leading-order correction ATc = Tc — Tq to the transi- 
tion temperature of a weakly interacting uniform Bose 
gas; see refs. 1 and 2 for an overview. However, the 
results differ substantially among different approaches, 
with ATc ranging from negative to positive values. Un- 
like the quantities at T = 0, perturbational calcula- 
tions are difficult for ATc. This situation also high- 
lights the absence of an established mean-field theory for 
Bose-Einstein condensates (BEC), corresponding to the 
Bardeen-Cooper-Schirieffer theory of superconductivity, 
with which one could estimate ATc easily. Thus, we still 
do not have a quantitative description of BEC at finite 
temperatures, particularly near Tc. 

As was reported by Hohenberg and Martin in 1965'^ 
and emphasized recently by Griffin,'* there are a couple 
of key words that can be used in formulating a mean- field 
theory for BEC: "conserving" and "gapless." Here, appli- 
cation of the conventional Wick-factorization procedure 
to the Hamiltonian of interacting Bose particles, known 
as the Hartree-Fock-Bogoliubov (HFB) theory,*' ^ fails to 
satisfy the Hugenholtz-Pines theorem^ giving rise to an 
unphysical energy gap in the excitation spectrum. Thus, 
it has become customary to introduce a further approxi- 
mation, now called the "Popov" approximation,* of com- 
pletely ignoring the anomalous quasiparticle pair corre- 
lation to recover a gapless excitation. However, when the 
same approximation is adopted for dynamics, it does not 
satisfy various conservation laws as required. Also, it is 
not clear whether it is permissible to neglect the pair 
correlation, which has played an essential role in the pio- 
neering perturbation theory of Bogoliubov;^ see also refs. 
8 and 9. Recently, Proukakis et al}^ have presented an 
improved gapless theory with the pair correlation. How- 
ever, it still does not satisfy the conservation laws. 

On the basis of these observations, we here formulate 
a new mean-field theory for BEC with the desired con- 
serving gapless character. As first shown by Baym*^ and 
discussed in detail by Hohenberg and Martin,^ the con- 
servation laws are well satisfied in derivative approx- 
imations" where the self-energy E is given as a deriva- 
tive of some functional <I> = $[G] with respect to Green's 



function G. This relationship between E and G holds 
exactly for Fermi systems, as shown by Luttinger and 
Ward.^^ Indeed, <I> was first introduced by Luttinger and 
Ward as a part of the exact thermodynamic functional 
Vt = Vl[G]. It was used later by Baym^* to give the cri- 
terion for obtaining dynamical equations with conser- 
vation laws. We hence ask the question here: Can we 
construct a Luttinger-Ward thermodynamic functional 
for BEC within a mean-field approximation that also 
satisfies the Hugcnholtz-Pines theorem? This is indeed 
possible as will be shown below. The predictions of the 
resulting mean-field theory for the uniform system will 
be presented later with an expression for ATc. We put 
h^k^ = \ throughout. 

As usual, we express the field operator ip{r) as a sum 
of the condensate wave function ^'(r) and the quasipar- 
ticle field We also adopt Beliaev's Green's function 
approach'^'^ and define our Matsubara Green's function 
in Nambu space as 



G(r,r',T) = -f3(T,0(r,T)0t(r')) 
where fa is the third Pauli matrix and 



0: 



(1) 



(2) 



The factor fa in eq. (I) is usually absent^'* but essential 
for the present purpose. Indeed, the eigenvalue problem 
for G~^ then becomes equivalent to the Bogoliubov-de 
Gennes equation for quasiparticles, as seen below in eqs. 
(6) and (f fa), so that the logarithmic term in eq. (3) can 
be written entirely with respect to the quasiparticles. 

Using the Fourier transform G(r, r', a;„) with ujn being 
the Matsubara frequency, we write down our Luttinger- 
Ward functional n = U.{'^ ,G) as 



j ■^*{r){Ho-^i)-^{r) dr 



2/3, 



^ Tr ln(i7o + E 



(3) 
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Here, Hq denotes the kinetic-energy operator plus the 
external potential, /U is the chemical potential, P = T~^, 
1 is the unit matrix, and Hq and l{LUn) are defined by 

, l{u)n) = ' 



Ho 



Ho 














g-iw„0+ 



(4) 

with 0+ an infinitesimal positive constant. The symbol Tr 
in eq. (3) also includes an integration over space variables 
with multiplications of Hq, 73, and 1 by the unit matrix 
6{i — r') implied. Finally, S is the irreducible self-energy 
matrix obtained from the functional $ = $(^', ^f*, G) by 

S(r, r', iJn) = -2/3 S^/dG{r', r, a;„) . (5) 

With this relation, fl is stationary with respect to varia- 
tion in G satisfying Dyson's equation: 



zw„l - Ho 



(6) 



A key quantity in eq. (3) is We choose it so that the 
conservation laws and the Hugenholtz-Pines theorem are 
simultaneously satisfied. Explicitly, it is given by 

$ = Jdrjdr' U{r-r') |i |4'(r)n*(r')P 



\^ir)\^-TrnG{r',r',cOn)lM 



+ -Trfs^{r)^\r')G{r',r,Wn)i{uJn) 



^2/322^ 



-Trf3G(r,r,a;„)l(w„) 



x-TrT3G{r',r',Un')i{oJn') 



,(7) 



whore U denotes the interaction potential and ^' and 
are defined in the same way as eq. (2). The first 
and second terms in the two square brackets of eq. (7) 
are the Hartree and Fock terms in Nambu space, re- 
spectively. The difference between this theory and the 
HFB theory lies in the Fock terms. Indeed, the HFB 
theory^' ^ can be reproduced from eq. (7) by replacing 
G and f3'5'(r)^'^(r') in the two Fock terms by fsG 
and *(r)'4'^(r'), respectively; the present functional was 
found by applying the reverse procedure. 

Now, the self-energy can be calculated explicitly with 
eq. (5), yielding the (1, 1) and (1,2) components as 

S(r,r') = d{r-r') J dr"U{r-r") [|*(r")|^ 

(r")0(r"))] 

+U{r-r') [*(r)**(r') + ((/>'f(r')(/>(r))] , (8a) 

A(r, r') = U{r-r') [*(r)*(r') + (^(r)(/)(r'))] , (8b) 

respectively. The (2, 1) and (2, 2) components are given 
by — A*(r,r') and — I]*(r,r'), respectively. We can also 
derive the equation for ^'(r) from Sfl/S^*{r) = 0. Noting 



5^1 /5G = 6, we only need to consider the explicit ^* 
dependences to obtain 

(Fo-M)*(r)+ j [S(r,r')*(r')-A(r,r')**(r')]rfr' = 0. 

(9) 

In the uniform case with no external potential where 
^ = constant, eq. (9) reduces to the Hugenholtz-Pines 

relation /i = Efc=o^Afe=o! desired. The expression for 
the particle number N = —dfl/djj, is obtained similarly 



as 



N : 



dr[|vl'(r)p + (0t(^)0(^))] 



(10) 



In order to calculate thermodynamic quantities, it is con- 
venient to diagonalize the Green's function of eq. (6). To 
this end, consider the following eigenvalue problem: 



H{r,r') 



w.(r') 
-vtir') 



dr' = 



-vtir) 



(11a) 



where H = Hq -|- S — /xfa , the subscript u specifies the 
eigenstate, and {u^jV^) should be normalized as 



[|u,(r)|2-|«,(r)|2] dr = l. 



(lib) 



One can show that is real when eq. (lib) can be 
satisfied.^ It also follows from eq. (9) that, without the 
condition (lib), eq. (11a) has a special solution E = 
for u = v = '^. We hence expect E^ > under the con- 
dition (lib) if the system is stable; the appearance of 
a negative eigenvalue evidences an instability of the as- 
sumed ^'(r). Finally, one can show by using the symme- 
try H = —TiH*Ti that the eigenstate of H corresponding 
to —Ei^ is given by {—v„, u*)"^. Now, we can provide ex- 
plicit expressions for ((/)'l'(r')(/)(r)) Gii(r, r', r = — 0+) 
and (0(r)0(r')) = -Gi2(r,r',O) as 

('/'t(r-')</'(^)>=E KM<(^')/.+^.(r-K(r')(l + ./.)] , 



{4>{r)4>{r')}: 



(12a) 

u,{r)v^{r')+v,ir)u,ir')] (0.5 + /,) , 

(12b) 

where f„ = (e^"/-^— 1) ^ is the Bose distribution function. 
In deriving eq. (12b), use has been made of the identity 
J2i, ^i^{''')v„{r') — J2i, v„{r)u^{r'), as can be proved from 
the completeness of the eigenfunctions from eq. (11). 

Equations (8)-(12) form a closed set of self-consistent 
equations satisfying both the Hugenholtz-Pines theorem 
and various conservation laws. Also, the pair correlation 
{^^) is adequately included in eqs. (9) and (11); neglect- 
ing this contribution yields the HFB-Popov theory. 

Using eqs. (8)-(12) and the procedure in §6 of ref. 14, 
it is possible to transform eq. (3) in equilibrium into an 
expression without Green's functions as 

Oeq = TEln(l-e-^''/^) -J2^- [ 

~\ j'^^j ^^'{^(^'^') [^{r'Wir) + {(l>^ {r)(t){r'))] 
-A*{r,r') [^{r'Mr) + (</.(r')0(r))] 1 . (13) 



J. Phys. Soc. Jpn. 



Full Paper 



Author Name 3 



The expression for the entropy can be obtained from 5* = 
—dQ/dT by differentiating eq. (3) in terms of the exphcit 
T dependences. It yields the well-known expression: 



[(1 + /.)ln(l + In /.] . (14) 



This completes the formulation of our mean-field theory. 

Now, let us apply eqs. (8)-(12) to a uniform system 
with no external potential and the interaction: 



U{r-r') = —d(r^r') 
m 



(15) 



where m is the particle mass and a is the s-wavc scatter- 
ing length. We first expand quantities dependent on a sin- 
gle argument r and two arguments (r.r') in \/-i/2g»fc 'r 
and y-^Q^'^ i'^-'^ )^ respectively, with V the volume. We 
then find from eq. (8) that the expansion coefficients Sfe 
and Afc are independent of fe; we hence put S 
and Afc A hereafter. Also, \E' = y/no {uq: conden- 
sate density) in eq. (9) so that the equation reduces to 
the Hugenholtz-Pines relation = E — A. Substituting 
these terms into eq. (11), we obtain Ek = ek{ek + 2A), 
Uk = ^/{£,k+Ek)/2Ek, and Vk = ^{^k-Ek)/2Ek with 
efc = k"^ /^.ra and = efe+A. Putting these terms back into 
eqs. (8) and (10), we obtain T, = 8TTna/m. with n = N/V 
and 

]_sr^ f ^k-Ek ^ 1 
y 2Ek Ek e^^/T-l 



no 



(16a) 



Ana 
m 



1 A 



1 + 



. (16b) 



The first term in the round brackets of eq. (16b) pro- 
duces an ultraviolet divergence. In accordance with eq. 
(15) where the low-energy scattering length a is used,^^ 
we now introduce an energy cutoff Cc to the integral 
in such a way that To <C Cc ^ l/raa?. Substituting 
To = 3.31n^/^/m, the condition is transformed into 1 <C 
Cc/Tq <^ 0.3(an^/^)~^, which can be well satisfied for 
an^/^^1. Other summations in cq. (16) are convergent, 
and eq. (16) is now simplified into 



(17a) 



1 - 



v'8ma' 



\ A, 
To 



/3 



no 



V27rAT 

37^52(A/T) 



= 2C(3/2f /^an^ 

n C(3/2)To^ 
where ^(3/2) is the Ricmann Q function, and 



, (17b) 



2^2x3/2 
V^C(3/2) 

\/27f 



'cosh2tsinhf , 1 
dt+- 



aX sinh 2t , 



C(3/2) 



y/x + 0{x) {x^Q), (18a) 



71" Jo 



sinh t 



aX sinh 2t . 



■ dt 



I - l.lQhs/x + 0{x) (a;^0). (18b) 



Note that both gi and §2 are given in powers of ^/x for 
X— >0, which causes anomalous behavior near in vari- 
ous thermodynamic quantities, as seen below. Equation 
(17) enables the complete determination of no and A, 
thereby fixing the thermodynamic equilibrium at a given 
temperature. The chemical potential is then obtained by 
the Hugcnholtz-Pines relation // = 8nna/m — A. Other 
quantities of interest are the entropy S of eq. (14) and 
the superfluid density Ps obtained by^ 



1 



1 



3mTV 



(e^fc/T _ 1)2 



(19) 



The transition temperature corresponds to the point 

where no = in cq. (17). Solving the rcsidting coupled 
equations for Tc and A(Tc) to the leading order in 5 = 
an}^^, we obtain 



To 



^ = 1 



A(T.) 



3C(3/2)V3 



87r 



1-^2.33(5, 



13.3(5^ 



(20a) 



(20b) 



To C(3/2)2/3 

Thus, the transition is first order with A(Tc) > 0; it is 
caused by the change in the low-energy quasiparticle dis- 
persion through Tf.. In contrast, no is continuous at in 
the present theory. Note the differences from the HFB- 
Popov theory where T^, = To and no is discontinuous at 
Tc. Our expression for ATc agrees with the analytic re- 
sult by Baym et al}^ as well as the numerical one by 
Holzmann and Krauth.^^ Quantities other than no are 
discontinuous at Tc as 

A/x(Tc) _ 2 A^(Tc) ^ 47r 

■ 3 N C(3/2)2/3' 



To 



6.63(5-' 



(20c) 



Ap,(Tc) 



47r 



= 1.16<5. 



(20d) 



nm 3C(3/2)V3 

Singular behaviors are also seen around T<Tc stemming 
from gi and (?2 in eq. (18). For example, 

no \/67r 



C(3/2)V3 



^1/2(1 _ 7^/2^)1 



/2 



-to -to 



S_ 
N 



To 



3/2 



5C(5/2) 3A 



(21a) 
(21b) 
(21c) 



2C(3/2) 2T 

It follows from eqs. (21b) and (21c) that the specific heat 

C = T{dS/dT) just below Tc is divergent as C/N - 
13.0(53/2(i-T/rcj-i/2_ 

Similar calculations at T = lead to the following ex- 
pressions in agreement with the Bogoliubov theory:^ 

«o(0) , 



30F 



^3/2 = 1 - 1.50^3/2 . 



m A(o) 



= 2C(3/2)^/^^ = 3.79(5. 



(22a) 



(22b) 



To To 

Figures 1-3 show the overall temperature dependences 
of C/N, no/n, ps/nm, and h/Tq as functions of T/Tc 
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0.0 0.2 0.4 0.6 0.8 1.0 

Fig. 1. Specific heat C/N as a function of T/Tc for an^/^ =0.05 
shown in comparison with the ideal gas result. 




0.0 0.2 0.4 0.6 

777 



Fig. 2. Normalized condensate density no/n and superfluid den- 
sity ps/nm as functions of T/Tc for an^^^ =0.05 shown in com- 
parison with the ideal gas result for no/n {=ps/nm). 



for a weakly interacting Bose gas with a-n}/^ — Q.Qb. For 
comparison, the corresponding results for the ideal gas 
(a = 0; Tc = Tq) are also plotted. The curves for av}/^ = 
0.05 satisfy the limiting behaviors derived above for T— s- 
Tc and 0, although Aps(7c) and 1— no(0)/n are too small 
on the present scale to be seen clearly. 

An important feature of the present mean-field theory 
is that the thermodynamic equilibrium is fixed by the 
coupled eq. (f 7) for tiq and A. Thus, the pair correlation 
{ipcj)) also plays an essential role for Tc- Now, one may 
ask; Arc the results for the transition compatible with the 
second-order transition in superfluid *He with Tc < Tq? It 
is expected that, as the interaction between particles is 
increased from zero, Tc initially increases towards a max- 
imum, decreasing eventually below Tq, as rationalized by 
Griiter et al}^ Also, the transition may change its char- 
acter during the course from first- to second-order. 

We finally comment on our result ATc = 2.3San^^^To in 
connection with preceding numerical estimates. Among 
them, the path-integral approach by Holzmann and 
Krauth^^ seems the most reliable, since they calculated 
ATc in terms of a correlation function in the noninter- 
acting Bose gas derived by an expansion from a = 0. 
They thereby obtained the result ATc^2.3an^/'^To com- 
patible with our analytical one. On the other hand, a 
different value ATc^l.San^^'^To has been obtained from 



Monte Carlo calculations by Arnold and Moore^ and by 

0.06, ■ , ■ , ' , ' , ' r- 



0.04 



0.02 



0.00 



an = 0.05 



ideal 



0.0 0.2 0.4 0.6 0.8 1.0 
T/T 



Fig. 3. Temperature dependence of the chemical potential fi for 
an^/^ =0.05 shown in comparison with the ideal gas result. 



Kashurnikov et al.^^ However, they both adopted as a 
starting point an effective action with a finite interac- 
tion, where the effects from the pair correlation {(fxj)) may 
not have been included appropriately. Finally, it is worth 
mentioning that the higher-order corrections to eq. (20a) 
depend on the cut-off Cc and are hence model dependent. 

In summary, we have derived new mean-field equa- 
tions for BEC as eqs. (8)-(12), which are appHcable to 
trapped atomic gases at finite temperatures. They have 
been applied to a uniform system to reveal its basic ther- 
modynamic features, particularly around Tc- 
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